import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
!pwd
/Users/tony/project_planckVsWMAP
hp.read_map("ffp10_newdust_total_353_full_map.fits") #planckTmap
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:368: UserWarning: If you are not specifying the input dtype and using the default np.float64 dtype of read_map(), please consider that it will change in a future version to None as to keep the same dtype of the input file: please explicitly set the dtype if it is important to you.
warnings.warn(
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:391: UserWarning: NSIDE = 2048
warnings.warn("NSIDE = {0:d}".format(nside))
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:400: UserWarning: ORDERING = RING in fits file
warnings.warn("ORDERING = {0:s} in fits file".format(ordering))
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:428: UserWarning: INDXSCHM = IMPLICIT
warnings.warn("INDXSCHM = {0:s}".format(schm))
array([ 2.50327867e-06, -4.45016194e-05, 3.31364572e-05, ...,
-4.48527280e-05, -1.91247091e-05, 2.42725946e-05])
map_sample = hp.read_map("ffp10_newdust_total_030_full_map.fits")
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:391: UserWarning: NSIDE = 1024
warnings.warn("NSIDE = {0:d}".format(nside))
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:400: UserWarning: ORDERING = NESTED in fits file
warnings.warn("ORDERING = {0:s} in fits file".format(ordering))
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:486: UserWarning: Ordering converted to RING
warnings.warn("Ordering converted to RING")
map_sample #353Hz map
array([-2.71337165e-04, -9.36984798e-05, -2.51075049e-04, ...,
-1.52736000e-04, -1.10948109e-04, -7.72219792e-05])
hp.mollview(
map_sample,
coord=["G", "E"],
title="Histogram equalized Ecliptic Sample 30Hz",
unit="mK",
norm="hist",
min=-0.01,
max=0.01,
)
hp.graticule()
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:920: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
newcm.set_over(newcm(1.0))
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:921: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
newcm.set_under(bgcolor)
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:922: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
newcm.set_bad(badcolor)
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:202: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
aximg = self.imshow(
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:541: UserWarning: 0.0 180.0 -180.0 180.0
warnings.warn(
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:656: UserWarning: The interval between parallels is 30 deg -0.00'.
warnings.warn(
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:664: UserWarning: The interval between meridians is 30 deg -0.00'.
warnings.warn(
Beam Maps
import astropy.io.fits as pf
import pylab # only to produce plots
FITSfile = 'BeamWf_HFI_R3.01/Bl_T_R3.01_fullsky_353x353.fits'
pf.info(FITSfile) # print list of extensions found in FITSfile
data, header = pf.getdata(FITSfile, header=True)#, 10, header=True) # read extension #10 (data and header)
#data, header = pf.getdata(FITSfile100, '100', header=True) # read extension having EXTNAME='ABC' (data and header)
print(header) # print header
print(data.names) # print column names
pylab.plot( data.field(0).flatten() ) # plot 1st column of binary table
Filename: BeamWf_HFI_R3.01/Bl_T_R3.01_fullsky_353x353.fits No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 19 () 1 WINDOW FUNCTION 1 TableHDU 26 4001R x 1C [E15.7] XTENSION= 'TABLE ' / ASCII table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 15 / length of dimension 1 NAXIS2 = 4001 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 1 / number of table fields TTYPE1 = 'TEMPERATURE' TFORM1 = 'E15.7 ' TBCOL1 = 1 EXTNAME = 'WINDOW FUNCTION' / extension name MAX-LPOL= 4000 / Maximum L multipole POLAR = F BCROSS = F ASYMCL = F COMMENT ---------------------------------------------------------------- COMMENT Beam Window Function B(l) COMMENT Compatible with Healpix (synfast, smoothing, ...) and PolSpice COMMENT To be squared before applying to power spectrum COMMENT C_map(l) = C_sky(l) * B(l)^2 COMMENT Adapted from COMMENT /redtruck/?????/m3space/quickbeam/quickpol/data/RD12rc3/beam_matrix_353xCOMMENT 353_l4000_s6_A000_cmbfast_0_rbIMO_rhIMO.npz COMMENT by ./qp2fits.py on 2016-06-15 COMMENT ---------------------------------------------------------------- END ['TEMPERATURE']
[<matplotlib.lines.Line2D at 0x102eca520>]
data
FITS_rec([(1.0,), (0.99999952,), (0.99999863,), ..., (0.068820141,),
(0.068731263,), (0.068642475,)],
dtype=(numpy.record, [('TEMPERATURE', 'S15')]))
#Intensity Map (Planck) - Foreground removed
intensityCMB100 = hp.read_map("COM_CMB_IQU-100-fgsub-sevem-field-Int_2048_R2.01_full.fits")
intensityCMB100
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:368: UserWarning: If you are not specifying the input dtype and using the default np.float64 dtype of read_map(), please consider that it will change in a future version to None as to keep the same dtype of the input file: please explicitly set the dtype if it is important to you.
warnings.warn(
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:391: UserWarning: NSIDE = 2048
warnings.warn("NSIDE = {0:d}".format(nside))
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:400: UserWarning: ORDERING = NESTED in fits file
warnings.warn("ORDERING = {0:s} in fits file".format(ordering))
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:426: UserWarning: No INDXSCHM keyword in header file : assume IMPLICIT
warnings.warn("No INDXSCHM keyword in header file : " "assume {}".format(schm))
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:428: UserWarning: INDXSCHM = IMPLICIT
warnings.warn("INDXSCHM = {0:s}".format(schm))
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/fitsfunc.py:486: UserWarning: Ordering converted to RING
warnings.warn("Ordering converted to RING")
array([-1.34818620e-04, -1.13565744e-04, -5.30425641e-05, ...,
9.37869918e-05, 3.33317512e-05, -2.04297194e-05])
hp.mollview(
intensityCMB100,
coord=["G", "E"],
title="Intensity CMB",
unit="mK",
norm="hist",
min=-0.01,
max=0.01,
)
hp.graticule()
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:920: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
newcm.set_over(newcm(1.0))
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:921: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
newcm.set_under(bgcolor)
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:922: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
newcm.set_bad(badcolor)
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:202: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
aximg = self.imshow(
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:541: UserWarning: 0.0 180.0 -180.0 180.0
warnings.warn(
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:656: UserWarning: The interval between parallels is 30 deg -0.00'.
warnings.warn(
/opt/anaconda3/envs/act_notebooks/lib/python3.9/site-packages/healpy/projaxes.py:664: UserWarning: The interval between meridians is 30 deg -0.00'.
warnings.warn(
intensityCMB217 = hp.read_map("COM_CMB_IQU-217-fgsub-sevem-field-Int_2048_R2.01_full.fits")
intensityCMB217
array([-1.45112805e-04, -9.79054021e-05, -8.30674253e-05, ...,
1.03153470e-04, 9.47897715e-05, 8.28015836e-05])
hp.mollview(
intensityCMB217,
coord=["G", "E"],
title="Intensity CMB",
unit="mK",
norm="hist",
min=-0.01,
max=0.01,
)
hp.graticule()
WMAP
!pwd
/Users/tony/project_planckVsWMAP
import healpy as hp
import numpy as np
import os
import astropy.units as u
import matplotlib.pyplot as plt
%matplotlib inline
.
Tool used
!curl "https://irsa.ipac.caltech.edu/data/Planck/release_3/all-sky-maps/maps/component-maps/cmb/COM_CMB_IQU-commander_2048_R3.00_full.fits" -o COM_CMB_IQU-commander_2048_R3.00_full.fits
filename = 'COM_CMB_IQU-commander_2048_R3.00_full.fits'
cmb_map = hp.read_map(filename)
Tool used
!curl "https://irsa.ipac.caltech.edu/data/Planck/release_3/ancillary-data/masks/COM_Mask_CMB-common-Mask-Int_2048_R3.00.fits" -o COM_Mask_CMB-common-Mask-Int_2048_R3.00.fits
#Visualising cmb_map
hp.mollview(cmb_map, min=-1e-3, max=1e-3, title="CMB only temperature map", unit="K")
Observation: There is residual galactic emission we should mask. The contamination just close to the galactic plane, (we could run anafast and specify a few degrees of gal_cut).
#masking with one of the planck masks(earlier loaded)
path = 'COM_Mask_CMB-common-Mask-Int_2048_R3.00.fits'
mask = hp.read_map(path)
map_masked = hp.ma(cmb_map)
map_masked.mask = np.logical_not(mask)
mask
array([1., 1., 1., ..., 1., 1., 1.])
map_masked
masked_array(data=[-0.00014912881306372583, -0.00011455055937403813,
-9.056342241819948e-05, ..., 0.00011307044769637287,
9.498728468315676e-05, 0.00010467269021319225],
mask=[False, False, False, ..., False, False, False],
fill_value=-1.6375e+30)
map_masked.mask
array([False, False, False, ..., False, False, False])
hp.mollview(map_masked, min=-1e-3, max=1e-3)
Now,we can load the binned TT CMB power spectrum:
Tool used
!curl "https://irsa.ipac.caltech.edu/data/Planck/release_3/ancillary-data/cosmoparams/COM_PowerSpect_CMB-TT-binned_R3.01.txt" -o COM_PowerSpect_CMB-TT-binned_R3.01.txt
!head -3 COM_PowerSpect_CMB-TT-binned_R3.01.txt #use of header info?
# l Dl -dDl +dDl BestFit 4.77112240e+01 1.47933552e+03 5.07654876e+01 5.07654876e+01 1.46111304e+03 7.64716065e+01 2.03496833e+03 5.47101576e+01 5.47101576e+01 2.06238073e+03
cmb_binned_spectrum = np.loadtxt('COM_PowerSpect_CMB-TT-binned_R3.01.txt')
cmb_binned_spectrum
array([[4.77112240e+01, 1.47933552e+03, 5.07654876e+01, 5.07654876e+01,
1.46111304e+03],
[7.64716065e+01, 2.03496833e+03, 5.47101576e+01, 5.47101576e+01,
2.06238073e+03],
[1.05917385e+02, 2.95539416e+03, 6.49766440e+01, 6.49766440e+01,
2.90452526e+03],
[1.35605348e+02, 3.86951392e+03, 7.69143744e+01, 7.69143744e+01,
3.90248963e+03],
[1.65405597e+02, 4.88946506e+03, 8.65856259e+01, 8.65856259e+01,
4.86135537e+03],
[1.95266870e+02, 5.46410945e+03, 9.05533459e+01, 9.05533459e+01,
5.53480715e+03],
[2.25164945e+02, 5.79343954e+03, 8.71348811e+01, 8.71348811e+01,
5.72693135e+03],
[2.55086908e+02, 5.37288375e+03, 7.69383697e+01, 7.69383697e+01,
5.37498113e+03],
[2.85025248e+02, 4.62767753e+03, 6.25199141e+01, 6.25199141e+01,
4.58832266e+03],
[3.14975304e+02, 3.60423851e+03, 4.71589439e+01, 4.71589439e+01,
3.59689801e+03],
[3.44934027e+02, 2.63120029e+03, 3.38762191e+01, 3.38762191e+01,
2.66745188e+03],
[3.74899344e+02, 2.03305943e+03, 2.48190562e+01, 2.48190562e+01,
2.01816864e+03],
[4.04869791e+02, 1.75336253e+03, 2.07109083e+01, 2.07109083e+01,
1.74834025e+03],
[4.34844309e+02, 1.78757901e+03, 2.06713695e+01, 2.06713695e+01,
1.82192489e+03],
[4.64822111e+02, 2.16204649e+03, 2.27523861e+01, 2.27523861e+01,
2.09823728e+03],
[4.94802601e+02, 2.42208480e+03, 2.49726104e+01, 2.49726104e+01,
2.39772998e+03],
[5.24785320e+02, 2.57348050e+03, 2.59184118e+01, 2.59184118e+01,
2.57144487e+03],
[5.54769905e+02, 2.54629768e+03, 2.50510538e+01, 2.50510538e+01,
2.55180645e+03],
[5.84756069e+02, 2.36064530e+03, 2.27209283e+01, 2.27209283e+01,
2.36436140e+03],
[6.14743583e+02, 2.09543505e+03, 1.98505488e+01, 1.98505488e+01,
2.10480607e+03],
[6.44732258e+02, 1.88467698e+03, 1.75025264e+01, 1.75025264e+01,
1.89096016e+03],
[6.74721939e+02, 1.81316288e+03, 1.64306044e+01, 1.64306044e+01,
1.81252216e+03],
[7.04712498e+02, 1.88319393e+03, 1.67990643e+01, 1.67990643e+01,
1.89672650e+03],
[7.34703827e+02, 2.09713072e+03, 1.81645476e+01, 1.81645476e+01,
2.10225612e+03],
[7.64695836e+02, 2.31873584e+03, 1.97363538e+01, 1.97363538e+01,
2.33880907e+03],
[7.94688447e+02, 2.46458129e+03, 2.07062220e+01, 2.07062220e+01,
2.50417440e+03],
[8.24681596e+02, 2.52191260e+03, 2.05197620e+01, 2.05197620e+01,
2.52277302e+03],
[8.54675226e+02, 2.39412077e+03, 1.90381031e+01, 1.90381031e+01,
2.37141912e+03],
[8.84669287e+02, 2.08333867e+03, 1.65790348e+01, 1.65790348e+01,
2.08335227e+03],
[9.14663738e+02, 1.74072841e+03, 1.37374917e+01, 1.37374917e+01,
1.73253226e+03],
[9.44658541e+02, 1.41867342e+03, 1.11420483e+01, 1.11420483e+01,
1.40502882e+03],
[9.74653663e+02, 1.17295424e+03, 9.28646475e+00, 9.28646475e+00,
1.17014623e+03],
[1.00464908e+03, 1.06240462e+03, 8.37987000e+00, 8.37987000e+00,
1.05887288e+03],
[1.03464476e+03, 1.04786241e+03, 8.30074773e+00, 8.30074773e+00,
1.05936420e+03],
[1.06464068e+03, 1.13235191e+03, 8.68575434e+00, 8.68575434e+00,
1.12669220e+03],
[1.09463683e+03, 1.21101851e+03, 9.11274984e+00, 9.11274984e+00,
1.20298244e+03],
[1.12463318e+03, 1.23187224e+03, 9.26050628e+00, 9.26050628e+00,
1.23858790e+03],
[1.15462972e+03, 1.20578050e+03, 8.98528122e+00, 8.98528122e+00,
1.20804315e+03],
[1.18462643e+03, 1.11792662e+03, 8.33171031e+00, 8.33171031e+00,
1.11502109e+03],
[1.21462331e+03, 9.68422243e+02, 7.54508665e+00, 7.54508665e+00,
9.86438951e+02],
[1.24462034e+03, 8.64397967e+02, 6.64322644e+00, 6.64322644e+00,
8.59469330e+02],
[1.27461751e+03, 7.61605936e+02, 5.98585735e+00, 5.98585735e+00,
7.66845300e+02],
[1.30461481e+03, 7.32339633e+02, 5.68666086e+00, 5.68666086e+00,
7.25319871e+02],
[1.33461223e+03, 7.37290337e+02, 5.71800485e+00, 5.71800485e+00,
7.31898515e+02],
[1.36460976e+03, 7.74925851e+02, 5.94089388e+00, 5.94089388e+00,
7.67047643e+02],
[1.39460740e+03, 8.06625346e+02, 6.17243469e+00, 6.17243469e+00,
8.03588527e+02],
[1.42460514e+03, 8.09057064e+02, 6.25632882e+00, 6.25632882e+00,
8.16784389e+02],
[1.45460297e+03, 7.77751089e+02, 6.11143430e+00, 6.11143430e+00,
7.92658706e+02],
[1.48460089e+03, 7.28943580e+02, 5.74544667e+00, 5.74544667e+00,
7.31230681e+02],
[1.51459890e+03, 6.48288645e+02, 5.24258624e+00, 5.24258624e+00,
6.44972815e+02],
[1.54459698e+03, 5.51286111e+02, 4.72393254e+00, 4.72393254e+00,
5.52838213e+02],
[1.57459513e+03, 4.76600995e+02, 4.30079600e+00, 4.30079600e+00,
4.73410287e+02],
[1.60459335e+03, 4.19541270e+02, 4.04442416e+00, 4.04442416e+00,
4.18888542e+02],
[1.63459164e+03, 3.95211496e+02, 3.96833591e+00, 3.96833591e+00,
3.92152600e+02],
[1.66458999e+03, 3.91657847e+02, 4.03423921e+00, 4.03423921e+00,
3.87313797e+02],
[1.69458839e+03, 3.92947519e+02, 4.17356282e+00, 4.17356282e+00,
3.93097715e+02],
[1.72458686e+03, 3.97735803e+02, 4.31665510e+00, 4.31665510e+00,
3.97391812e+02],
[1.75458537e+03, 3.83011743e+02, 4.41586492e+00, 4.41586492e+00,
3.91427360e+02],
[1.78458394e+03, 3.75391143e+02, 4.45733427e+00, 4.45733427e+00,
3.71983403e+02],
[1.81458255e+03, 3.40139516e+02, 4.45892926e+00, 4.45892926e+00,
3.41585903e+02],
[1.84458121e+03, 3.07539268e+02, 4.45749277e+00, 4.45749277e+00,
3.06532670e+02],
[1.87457991e+03, 2.73859536e+02, 4.49350451e+00, 4.49350451e+00,
2.74153525e+02],
[1.90457865e+03, 2.49846901e+02, 4.59652807e+00, 4.59652807e+00,
2.50006349e+02],
[1.93457743e+03, 2.44246300e+02, 4.77758093e+00, 4.77758093e+00,
2.36236579e+02],
[1.96457625e+03, 2.31913142e+02, 5.02866666e+00, 5.02866666e+00,
2.31349665e+02],
[1.99457510e+03, 2.34658998e+02, 5.42839879e+00, 5.42839879e+00,
2.31209089e+02],
[2.02457399e+03, 2.31744942e+02, 5.86880339e+00, 5.86880339e+00,
2.30818750e+02],
[2.05457291e+03, 2.19660112e+02, 6.20024969e+00, 6.20024969e+00,
2.26094532e+02],
[2.08457186e+03, 2.14252650e+02, 6.53323083e+00, 6.53323083e+00,
2.15102553e+02],
[2.11457084e+03, 1.96716713e+02, 6.87393629e+00, 6.87393629e+00,
1.98340298e+02],
[2.14456985e+03, 1.76431323e+02, 7.23789713e+00, 7.23789713e+00,
1.78237039e+02],
[2.17456889e+03, 1.52227062e+02, 7.64417922e+00, 7.64417922e+00,
1.58049842e+02],
[2.20456795e+03, 1.38671109e+02, 8.10965997e+00, 8.10965997e+00,
1.40733305e+02],
[2.23456704e+03, 1.21459219e+02, 8.64541703e+00, 8.64541703e+00,
1.28004092e+02],
[2.26456615e+03, 1.16155153e+02, 9.25556670e+00, 9.25556670e+00,
1.20051892e+02],
[2.29456529e+03, 1.17768892e+02, 9.93861339e+00, 9.93861339e+00,
1.15705770e+02],
[2.32456444e+03, 1.16129438e+02, 1.06905048e+01, 1.06905048e+01,
1.13100669e+02],
[2.35456362e+03, 1.09714703e+02, 1.15085752e+01, 1.15085752e+01,
1.10345705e+02],
[2.38456282e+03, 1.20361413e+02, 1.23926722e+01, 1.23926722e+01,
1.06139560e+02],
[2.41456204e+03, 9.04515698e+01, 1.33433129e+01, 1.33433129e+01,
1.00105444e+02],
[2.44456128e+03, 1.03080972e+02, 1.43819739e+01, 1.43819739e+01,
9.26666904e+01],
[2.47456054e+03, 7.36743738e+01, 1.55113561e+01, 1.55113561e+01,
8.48106882e+01],
[2.49902400e+03, 6.13290585e+01, 2.09398269e+01, 2.09398269e+01,
7.87736295e+01]])
lmax = 3000
test_cls_meas_frommap = hp.anafast(map_masked, lmax=lmax, use_pixel_weights=True) #used "use_pixel_weights=True" arg to have a more precise power spectrum
ll = np.arange(lmax+1)
ll
array([ 0, 1, 2, ..., 2998, 2999, 3000])
sky_fraction = len(map_masked.compressed()) / len(map_masked)
print(f"The map covers {sky_fraction:.1%} of the sky")
The map covers 77.9% of the sky
plt.style.use("seaborn-poster") #read more on style used
k2muK = 1e6
Power spectra are generally plotted as $D_\ell$ which is defined as $\dfrac{\ell(\ell+1)}{2 \pi}C_\ell$, so we need to apply that factor to the $C_\ell$ calculated from the map.
plt.plot(cmb_binned_spectrum[:,0], cmb_binned_spectrum[:,1], '--', alpha=1, label='Planck 2018 PS release')
plt.plot(ll, ll*(ll+1.)*test_cls_meas_frommap*k2muK**2/2./np.pi / sky_fraction, '--', alpha=0.6, label='Planck 2018 PS from Data Map')
plt.xlabel(r'$\ell$')
plt.ylabel(r'$D_\ell~[\mu K^2]$')
plt.grid()
plt.legend(loc='best')
<matplotlib.legend.Legend at 0x10ddcb160>
Reading the documentation of the Planck release we see that the output has a resolution of 5 arcminutes. Therefore as a first order correction of the beam, we can divide the power spectrum by the square of the beam window function.
w_ell = hp.gauss_beam((5*u.arcmin).to_value(u.radian), lmax=lmax)
plt.plot(cmb_binned_spectrum[:,0], cmb_binned_spectrum[:,1], '--', alpha=1, label='Planck 2018 PS release')
plt.plot(ll, ll*(ll+1.)*test_cls_meas_frommap*k2muK**2/2./np.pi / sky_fraction / w_ell**2,
alpha=0.6, label='Planck 2018 PS from Data Map (beam corrected)')
plt.xlabel(r'$\ell$')
plt.ylabel(r'$D_\ell~[\mu K^2]$')
plt.grid()
plt.legend(loc='best');
Intermediate(Getting Clean Map used for analysis)
from pixell import enmap, utils
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
from astropy.convolution import Gaussian2DKernel
from scipy.signal import convolve as scipy_convolve
from astropy.convolution import convolve
# Load the data
filename = get_pkg_data_filename('ffp10_newdust_total_030_full_map.fits')
hdu = fits.open(filename)[0]
hdu
<astropy.io.fits.hdu.image.PrimaryHDU at 0x142386670>
from astropy.convolution import Gaussian1DKernel
#the convolution module also includes built-in kernels that can be imported as
gauss = Gaussian1DKernel(stddev=2)
gauss.array
array([6.69151129e-05, 4.36341348e-04, 2.21592421e-03, 8.76415025e-03,
2.69954833e-02, 6.47587978e-02, 1.20985362e-01, 1.76032663e-01,
1.99471140e-01, 1.76032663e-01, 1.20985362e-01, 6.47587978e-02,
2.69954833e-02, 8.76415025e-03, 2.21592421e-03, 4.36341348e-04,
6.69151129e-05])
The kernel can then be used directly when calling convolve()
NSIDE = 1024
print("Approximate resolution at NSIDE {} is {:.2} deg".format(
NSIDE, hp.nside2resol(NSIDE, arcmin=True) / 60
) )
Approximate resolution at NSIDE 1024 is 0.057 deg
NPIX = hp.nside2npix(NSIDE)
print(NPIX)
12582912
getting alms of map(sample) used from Planck
hp.sphtfunc.map2alm(map_sample, lmax=1300, mmax=None, iter=3, pol=True, use_weights=False, datapath=None, gal_cut=0, use_pixel_weights=False)
array([ 1.16211765e-03+0.00000000e+00j, -1.29597927e-05+0.00000000e+00j,
-1.32805698e-03+0.00000000e+00j, ...,
2.49503412e-08-4.73206538e-08j, -1.51505776e-08-3.88394977e-08j,
-2.92800850e-08-3.33140654e-08j])
#pol=false
alm = hp.sphtfunc.map2alm(map_sample, lmax=1300, mmax=None, iter=3, pol=False, use_weights=False, datapath=None, gal_cut=0, use_pixel_weights=False)
alm
array([ 1.16211765e-03+0.00000000e+00j, -1.29597927e-05+0.00000000e+00j,
-1.32805698e-03+0.00000000e+00j, ...,
2.49503412e-08-4.73206538e-08j, -1.51505776e-08-3.88394977e-08j,
-2.92800850e-08-3.33140654e-08j])
" Multiply alm by a function of l. The function is assumed to be zero where not defined.
Parameters almarray The alm to multiply
flarray The function (at l=0..fl.size-1) by which alm must be multiplied.
mmaxNone or int, optional The maximum m defining the alm layout. Default: lmax.
inplacebool, optional If True, modify the given alm, otherwise make a copy before multiplying.
Returns almarray The modified alm, either a new array or a reference to input alm, if inplace is True. "
np.shape(alm)
(846951,)
alm
array([ 1.16211765e-03+0.00000000e+00j, -1.29597927e-05+0.00000000e+00j,
-1.32805698e-03+0.00000000e+00j, ...,
2.49503412e-08-4.73206538e-08j, -1.51505776e-08-3.88394977e-08j,
-2.92800850e-08-3.33140654e-08j])
#fl = w_ell #see power spectra step, deriv of beam function used as fl (function of l)
#hp.sphtfunc.almxfl(alm , fl, mmax=None, inplace=True)
alm
array([ 1.16211765e-03+0.00000000e+00j, -1.29597927e-05+0.00000000e+00j,
-1.32805698e-03+0.00000000e+00j, ...,
2.49503412e-08-4.73206538e-08j, -1.51505776e-08-3.88394977e-08j,
-2.92800850e-08-3.33140654e-08j])
applying a spherical harmonic transform back to map space (hp.alm2map), using Nside = 1024.
#hp.sphtfunc.alm2map(alm, nside=1024, lmax=lmax, mmax=None, pixwin=False, fwhm=0.0, sigma=None, pol=True, inplace=False, verbose=True)
#mmax = 1300
#checking size of alm
(mmax * (2 * lmax + 1 - mmax)) / (2 + lmax + 1)
2035.064935064935
#resize alm
#alm_resize= np.resize(alm,2035)
#alm_resize
#pol=false polarisatn
#hp.sphtfunc.alm2map(alm,nside=16, lmax=lmax, mmax=None, pixwin=False, fwhm=0.0, sigma=None, pol=False, inplace=False, verbose=True)
new_map=hp.sphtfunc.alm2map(alm,nside=1024) #, mmax=None, pixwin=False, fwhm=0.0, sigma=None, pol=False, inplace=False, verbose=True)
new_map
array([-0.00018566, -0.00016411, -0.00019511, ..., -0.00013042,
-0.0001177 , -0.0001039 ])
np.shape(new_map)
(12582912,)
np.shape(alm)
(846951,)
np.shape(map_sample)
(12582912,)
Alternative approach for alm2map
hp.sphtfunc.alm2cl(alm, alms2=None, lmax=None, mmax=None, lmax_out=None, nspec=None)
#Computes (cross-)spectra from alm(s). If alm2 is given, cross-spectra between alm and alm2 are computed. If alm (and alm2 if provided) contains n alm, then n(n+1)/2 auto and cross-spectra are returned
array([1.35051743e-06, 2.53453754e-07, 3.84581273e-07, ...,
3.14515626e-15, 3.22686188e-15, 3.20514566e-15])
cl = hp.sphtfunc.alm2cl(alm, alms2=None, lmax=None, mmax=None, lmax_out=None, nspec=None)
cl
array([1.35051743e-06, 2.53453754e-07, 3.84581273e-07, ...,
3.14515626e-15, 3.22686188e-15, 3.20514566e-15])
alm2map = hp.sphtfunc.synfast(cl, nside=1024, lmax=1300, mmax=None, alm=False, pol=True, pixwin=False, fwhm=0.0, sigma=None, new=False, verbose=True)
alm2map
array([0.00039491, 0.00053628, 0.00071575, ..., 0.0001484 , 0.00048541,
0.00060199])
Creates a map(s) from cl(s). Parameters cls [array or tuple of array] A cl or a list of cl (either 4 or 6, see synalm()) nside [int, scalar] The nside of the output map(s) lmax [int, scalar, optional] Maximum l for alm. Default: min of 3*nside-1 or length of the cls - 1 mmax [int, scalar, optional] Maximum m for alm. alm [bool, scalar, optional] If True, return also alm(s). Default: False. pol [bool, optional] If True, assumes input cls are TEB and correlation. Output will be TQU maps. (input must be 1, 4 or 6 cl’s) If False, fields are assumed to be described by spin 0 spherical harmonics. (input can be any number of cl’s) If there is only one input cl, it has no effect. Default: True. pixwin [bool, scalar, optional] If True, convolve the alm by the pixel window function. Default: False. fwhm [float, scalar, optional] The fwhm of the Gaussian used to smooth the map (applied on alm) [in radians] sigma [float, scalar, optional] The sigma of the Gaussian used to smooth the map (applied on alm) [in radians] Returns maps [array or tuple of arrays] The output map (possibly list of maps if polarized input). or, if alm is True, a tuple of (map,alm) (alm possibly a list of alm if polarized input)
#estimating covariance matrix
#np.cov(alm2map, y=None, rowvar=True, bias=False, ddof=None, fweights=None, aweights=None, dtype=None) NOT helpful
alm2map
array([0.00039491, 0.00053628, 0.00071575, ..., 0.0001484 , 0.00048541,
0.00060199])
hp.mollview(
new_map,
coord=["G"], #galactic only
title="Visualised sample",
unit="mK",
norm="hist",
min=-0.01,
max=0.01,
)
hp.graticule()
hp.mollview(
map_sample,
coord=["G"],
title="Histogram equalized Galactic Sample 30Hz",
unit="mK",
norm="hist",
min=-0.01,
max=0.01,
)
hp.graticule()
hp.mollview(
map_sample - new_map,
coord=["G"],
title="Map Difference Galactic Sample 30Hz",
unit="mK",
norm="hist",
#min=-0.0,
#max=0.01,
)
hp.graticule()